plot(stdres(bestmodel.lmer) ~ d$FoodTreatment)
stdres(bestmodel.lmer)
resid(bestmodel.lmer)
?residuals.merMod
plot(resid(bestmodel.lmer, scaled=TRUE) ~ d$FoodTreatment)
10^0.168
10^-0.168
confint(bestmodel.lmer, oldNames = FALSE)
10^-0.207196073
10^-0.12942820
effects(bestmodel.lmer)
allEffects(bestmodel.lmer)
?allEffects
effect("FoodTreatment", bestmodel.lmer)
10^effect("FoodTreatment", bestmodel.lmer)
effect("FoodTreatment", bestmodel.lmer)
10^0.4115296
10^0.2431531
1.750464/2.579465
effect("ArrivalTime.centered", bestmodel.lmer)
?effect
allEffects(bestmodel.lmer)
??ggeffects
d <- read.delim(file = "Data/owls.txt")
d$Nest <- factor(d$Nest)
d$FoodTreatment <- factor(d$FoodTreatment)
d$SexParent <- factor(d$SexParent)
summary(d)
d <- read.delim(file = "Data/owls.txt")
d$Nest <- factor(d$Nest)
d$FoodTreatment <- factor(d$FoodTreatment)
d$SexParent <- factor(d$SexParent)
d$ArrivalTime.mean     <- ave(d$ArrivalTime, d$Nest)
d$ArrivalTime.centered <- d$ArrivalTime - d$ArrivalTime.mean
d$NegPerChick.log10    <- log10(d$NegPerChick + 1)
ggplot(data = d, aes(x = ArrivalTime.centered, y = NegPerChick.log10)) +
facet_wrap(~Nest) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Arrival time (centered for each nest)",
y = expression(log[10]("Negotiations per chick"))) +
theme_few()
library(ggplot2)
library(ggthemes)
ggplot(data = d, aes(x = ArrivalTime.centered, y = NegPerChick.log10)) +
facet_wrap(~Nest) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Arrival time (centered for each nest)",
y = expression(log[10]("Negotiations per chick"))) +
theme_few()
library(lme4)
randomslope.lmer <- lmer(
NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = d)
residplot(randomslope.lmer)
library(predictmeans)
randomslope.lmer <- lmer(
NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered | Nest), data = d)
residplot(randomslope.lmer)
randomintercept.lmer <- lmer(
NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(1 | Nest), data = d)
fixedeffects.lm <- lm(
NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3,
data = d)
randomslope.lmer.ml <- lmer(
NegPerChick.log10 ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered)^3 +
(ArrivalTime.centered || Nest), data = d, REML = FALSE, na.action = "na.fail")
library(MuMIn)
(dd <- subset(dredge(randomslope.lmer.ml), delta < 4))
plot(dd)
bestmodel.lmer <- lmer(
NegPerChick.log10 ~ FoodTreatment * ArrivalTime.centered + SexParent + (ArrivalTime.centered | Nest),
data = d)
S(bestmodel.lmer)
library(car)
S(bestmodel.lmer)
randomslope.glmer <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize))^3 +
(ArrivalTime.centered | Nest), data = d, family = "poisson")
randomslope.glmer <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize)))^3 +
(ArrivalTime.centered | Nest), data = d, family = "poisson")
randomslope.glmer <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize)))^2 +
(ArrivalTime.centered | Nest), data = d, family = "poisson")
randomslope.glmer <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize)))^3 +
(1 | Nest), data = d, family = "poisson")
randomslope.glmer <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize)))^3 +
(ArrivalTime.centered | Nest), data = d, family = "poisson")
S(g1)
S(randomslope.glmer )
randomslope.glmer2 <- glmer(
SiblingNegotiation ~ (SexParent + FoodTreatment + ArrivalTime.mean + ArrivalTime.centered + offset(log(BroodSize)))^3 +
(1 | Nest), data = d, family = "poisson")
REMLBoundaryTest <- function(fullerModel, reducedModel, df) {
D <- -2 * c((logLik(reducedModel, REML = TRUE) - logLik(fullerModel, REML = TRUE)))
p.value <- (pchisq(D, df, lower.tail = FALSE) + pchisq(D, df - 1, lower.tail = FALSE))/2
output <- c(D, df, p.value)
names(output) <- c("test statistic", "df", "p-value")
return(output)
}
REMLBoundaryTest(randomslope.glmer1, randomslope.glmer2, 2)
library(ggeffects)
?ggpredict
# Inspect the data
d <- read.csv(file = "Data/orthodont.csv")
d$Subject <- factor(d$Subject)
d$Sex <- factor(d$Sex)
summary(d)
ggpairs(d) +
theme_few()
library(car)
library(effects)
library(GGally)
library(ggplot2)
library(ggthemes)
library(lme4)
library(MuMIn)
library(predictmeans)
ggpairs(d) +
theme_few()
ggpairs(d[,-3]) +
theme_few()
ggplot(data = d, aes(x = age, y = distance, colour = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Arrival time (centered for each nest)",
y = expression(log[10]("Negotiations per chick"))) +
theme_few()
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Arrival time (centered for each nest)",
y = expression(log[10]("Negotiations per chick"))) +
theme_few()
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Arrival time (centered for each nest)",
y = expression(log[10]("Negotiations per chick"))) +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Age (years)", y = "Distance (mm)") +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
randomslope.lmer <- lmer(distance ~ Sex + (age | Subject), data = d)
residplot(randomslope.lmer)
ggpairs(d[,-3]) +
theme_few()
randomslope.lmer <- lmer(distance ~ Sex + age + (age | Subject), data = d)
residplot(randomslope.lmer)
randomslope.lmer <- lmer(distance ~ Sex * age + (age | Subject), data = d)
residplot(randomslope.lmer)
randomintercept.lmer <- lmer(distance ~ Sex * age + (1 | Subject), data = d)
fixedeffects.lm <- lm(distance ~ Sex * age, data = d)
REMLBoundaryTest(randomintercept.lmer, fixedeffects.lm, df = 1)
# Define REMLBoundaryTest() function for selecting random effects
REMLBoundaryTest <- function(fullerModel, reducedModel, df) {
D <- -2 * c((logLik(reducedModel, REML = TRUE) - logLik(fullerModel, REML = TRUE)))
p.value <- (pchisq(D, df, lower.tail = FALSE) + pchisq(D, df - 1, lower.tail = FALSE))/2
output <- c(D, df, p.value)
names(output) <- c("test statistic", "df", "p-value")
return(output)
}
REMLBoundaryTest(randomintercept.lmer, fixedeffects.lm, df = 1)
REMLBoundaryTest(randomslope.lmer, randomintercept.lmer, df = 2)
randomintercept.lmer.ml <- lmer(distance ~ Sex * age + (1 | Subject), data = d, na.action = "na.omit", REML = FALSE)
(dd <- subset(dredge(randomslope.lmer.ml), delta < 4))
randomintercept.lmer.ml <- lmer(distance ~ Sex * age + (1 | Subject), data = d, REML = FALSE, na.action = "na.omit")
(dd <- subset(dredge(randomslope.lmer.ml), delta < 4))
randomintercept.lmer.ml <- lmer(distance ~ Sex * age + (1 | Subject), data = d, REML = FALSE, na.action = "na.omit")
(dd <- subset(dredge(randomintercept.lmer.ml), delta < 4))
randomintercept.lmer.ml <- lmer(distance ~ Sex * age + (1 | Subject), data = d, REML = FALSE, na.action = "na.fail")
(dd <- subset(dredge(randomintercept.lmer.ml), delta < 4))
plot(dd)
residplot(randomintercept.lmer.ml)
# Re-fit the reduced model (using REML)
bestmodel.lmer <- lmer(distance ~ Sex * age + (1 | Subject), data = d)
# Check assumptions of the reduced model
residplot(bestmodel.lmer)
# Interpret the reduced model
S(bestmodel.lmer)
confint(bestmodel.lmer, oldNames = FALSE)
ranef(bestmodel.lmer)
r.squaredGLMM(bestmodel.lmer)
plot(allEffects(bestmodel.lmer),
main = NULL,
xlab = "Age (years)",
ylab = "Distance (mm)")
ggpairs(d[,-3]) +
theme_few()
ggpairs(d[,-3]) +
theme_few()
# Plot distance as a function of age with a panel for each subject, coloured by sex
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Age (years)", y = "Distance (mm)") +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
residplot(randomslope.lmer)
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Age (years)", y = "Distance (mm)") +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
ggpairs(d[,-3]) +
theme_few()
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Age (years)", y = "Distance (mm)") +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
randomslope.lmer <- lmer(distance ~ Sex * age + (age | Subject), data = d)
residplot(randomslope.lmer)
randomintercept.lmer <- lmer(distance ~ Sex * age + (1 | Subject), data = d)
fixedeffects.lm <- lm(distance ~ Sex * age, data = d)REMLBoundaryTest(randomintercept.lmer, fixedeffects.lm, df = 1)
# Test the significance of the random intercept
REMLBoundaryTest(randomintercept.lmer, fixedeffects.lm, df = 1)
REMLBoundaryTest(randomslope.lmer, randomintercept.lmer, df = 2)
# Model comparison for fixed effects (using ML)
# We will now re-fit the random-intercept model using ML and dredge the fixed effects
randomintercept.lmer.ml <- lmer(distance ~ Sex * age + (1 | Subject), data = d,
REML = FALSE, na.action = "na.fail")
(dd <- subset(dredge(randomintercept.lmer.ml), delta < 4))
plot(dd)
residplot(bestmodel.lmer)
S(bestmodel.lmer)
confint(bestmodel.lmer, oldNames = FALSE)
ranef(bestmodel.lmer)
r.squaredGLMM(bestmodel.lmer)
plot(allEffects(bestmodel.lmer),
main = NULL,
xlab = "Age (years)",
ylab = "Distance (mm)")
library(emmeans)
allEffects(bestmodel.lmer)
effect("age", bestmodel.lmer)
effect(Sex, bestmodel.lmer)
effect("Sex", bestmodel.lmer)
Anova(bestmodel.lmer, test = "F", type =3)
1.816^2
1.386^2
3.297856 + 1.920996
1.816^2 / (1.816^2 + 1.386^2)
S(bestmodel.lmer)
1.816^2 / (1.816^2 + 1.386^2)
1.816^2 / (1.816^2 + 1.386^2)
d <- read.csv(file = "Data/orthodont.csv")
d$Subject <- factor(d$Subject)
d$FoodTreatment <- factor(d$FoodTreatment)
d <- read.csv(file = "Data/orthodont.csv")
d$Subject <- factor(d$Subject)
d$FoodTreatment <- factor(d$FoodTreatment)
d$Sex <- factor(d$Sex)
summary(d)
# Center age
d$age.mean     <- ave(d$age, d$Subject)
d$age.centered <- d$age - d$age.mean
ggpairs(d[,-(2:3)]) +
theme_few()
# Tutorial 10 Answers
# Load packages
library(car)
library(effects)
library(GGally)
library(ggplot2)
library(ggthemes)
library(lme4)
library(MuMIn)
library(predictmeans)
# Define REMLBoundaryTest() function for selecting random effects
REMLBoundaryTest <- function(fullerModel, reducedModel, df) {
D <- -2 * c((logLik(reducedModel, REML = TRUE) - logLik(fullerModel, REML = TRUE)))
p.value <- (pchisq(D, df, lower.tail = FALSE) + pchisq(D, df - 1, lower.tail = FALSE))/2
output <- c(D, df, p.value)
names(output) <- c("test statistic", "df", "p-value")
return(output)
}
ggpairs(d[,-(2:3)]) +
theme_few()
randomslope.lmer <- lmer(distance ~ Sex * age.centered + (age.centered | Subject), data = d)
S(randomslope.lmer)
randomslope.lmer <- lmer(distance ~ Sex * age + (age | Subject), data = d)
S(randomslope.lmer)
randomslope.lmer <- lmer(distance ~ Sex * age.centered + (age.centered | Subject), data = d)
randomintercept.lmer <- lmer(distance ~ Sex * age.centered + (1 | Subject), data = d)
fixedeffects.lm <- lm(distance ~ Sex * age.centered, data = d)
REMLBoundaryTest(randomintercept.lmer, fixedeffects.lm, df = 1)
REMLBoundaryTest(randomslope.lmer, randomintercept.lmer, df = 2)
randomintercept.lmer.ml <- lmer(distance ~ Sex * age.centered + (1 | Subject), data = d,
REML = FALSE, na.action = "na.fail")
(dd <- subset(dredge(randomintercept.lmer.ml), delta < 4))
bestmodel.lmer <- lmer(distance ~ Sex * age.centered + (1 | Subject), data = d)
residplot(bestmodel.lmer)
S(bestmodel.lmer)
confint(bestmodel.lmer, oldNames = FALSE)
ranef(bestmodel.lmer)
r.squaredGLMM(bestmodel.lmer)
plot(allEffects(bestmodel.lmer),
main = NULL,
xlab = "Age (years)",
ylab = "Distance (mm)")
Anova(bestmodel.lmer, test = "F", type = 3)
bestmodel.lmer <- lmer(distance ~ Sex * age + (1 | Subject), data = d)
Anova(bestmodel.lmer, test = "F", type = 3)
S(bestmodel.lmer)
bestmodel.lmer <- lmer(distance ~ Sex * age.centered + (1 | Subject), data = d)
S(bestmodel.lmer)
S(bestmodel.lmer)
confint(bestmodel.lmer, oldNames = FALSE)
4.795e-01+3.048e-01
0.2964574+0.0669912
library(emmeans)
emtrends(model, ~ Sex, var = "age.centered")
library(emmeans)
emtrends(bestmodel.lmer, ~ Sex, var = "age.centered")
emmeans(bestmodel.lmer)
emmeans(bestmodel.lmer, ~Sex)
residuals.mermode
library(lme4)
?residuals
??resid
library(predictmeans)
?residplot
log10(1)
?residplot
# Tutorial 10 Answers
# Load packages
library(car)
library(effects)
library(GGally)
library(ggplot2)
library(ggthemes)
library(lme4)
library(MuMIn)
library(predictmeans)
output <- c(D, df, p.value)
names(output) <- c("test statistic", "df", "p-value")
# Define REMLBoundaryTest() function for selecting random effects
REMLBoundaryTest <- function(fullerModel, reducedModel, df) {
D <- -2 * c((logLik(reducedModel, REML = TRUE) - logLik(fullerModel, REML = TRUE)))
p.value <- (pchisq(D, df, lower.tail = FALSE) + pchisq(D, df - 1, lower.tail = FALSE))/2
output <- c(D, df, p.value)
names(output) <- c("test statistic", "df", "p-value")
return(output)
}
# Load the data, convert to factors, and inspect the data
d <- read.csv(file = "Data/orthodont.csv")
d$Subject <- factor(d$Subject)
d$Sex <- factor(d$Sex)
summary(d)
# Center age
d$age.mean     <- ave(d$age, d$Subject)
d$age.centered <- d$age - d$age.mean
# Plot the data as a scatterplot matrix (excluding the Subject column)
ggpairs(d[,-(2:3)]) +
theme_few()
# Plot distance as a function of age with a panel for each subject, coloured by sex
ggplot(data = d, aes(x = age, y = distance, colour = Sex, fill = Sex)) +
facet_wrap(~Subject) +
geom_point() +
geom_smooth(method = "lm") +
labs(x = "Age (years)", y = "Distance (mm)") +
scale_color_colorblind() +
scale_fill_colorblind() +
theme_few()
# Build the full model
# We will build a random-slope model including fixed effects of sex, age, and their interaction
# and random effects comprising a random slope and intercept for each subject
randomslope.lmer <- lmer(distance ~ Sex * age.centered + (age.centered | Subject), data = d)
# Check assumptions of the full model
residplot(randomslope.lmer, group = "Sex")
residplot(randomslope.lmer, group = "Sex", slope = TRUE)
randomslope.resid <- residuals(randomslope.lmer, scaled = TRUE)
plot(randomslope.resid ~ Sex, data = d)
plot(randomslope.resid ~ age.centered, data = d)
plot(randomslope.resid ~ Subject, data = d)
library(performance)
?binned
??binned
d <- read.csv("Data/storks.csv")
summary(d)
stork.glm <- glm(surv ~ cort, family = binomial, data = d, na.action = "na.fail")
binned_residuals(stork.glm)
residuals(stork.glm)
resid <- binned_residuals(stork.glm)
plot(resid)
plot(resid)
?see
check_model(stork.glm)
dredge(stork.glm)
library(MuMIn)
dredge(stork.glm)
d <- read.csv("Data/sheep.csv")
summary(d)
dredge(sheep.glm)
sheep.glm <- glm(fitness ~ body.size, family = poisson, data = d, na.action = "na.fail")
dredge(sheep.glm)
check_model(sheep.glm)
d <- read.delim(file = "Data/owls.txt")
d$Nest <- factor(d$Nest)
d$FoodTreatment <- factor(d$FoodTreatment)
d$SexParent <- factor(d$SexParent)
summary(d)
bestmodel.lmer <- lmer(
NegPerChick.log10 ~ FoodTreatment * ArrivalTime.centered + (ArrivalTime.centered | Nest),
data = d)
library(lme4)
bestmodel.lmer <- lmer(
NegPerChick.log10 ~ FoodTreatment * ArrivalTime.centered + (ArrivalTime.centered | Nest),
data = d)
d$ArrivalTime.mean     <- ave(d$ArrivalTime, d$Nest)
d$ArrivalTime.centered <- d$ArrivalTime - d$ArrivalTime.mean
d$NegPerChick.log10    <- log10(d$NegPerChick + 1)
bestmodel.lmer <- lmer(
NegPerChick.log10 ~ FoodTreatment * ArrivalTime.centered + (ArrivalTime.centered | Nest),
data = d)
check_model(bestmodel.lmer)
setwd("~/Library/CloudStorage/OneDrive-UniversityofNewBrunswick/Documents/Teaching/UNB Courses/2024-2025/Winter 2025/Experimental Design/Tutorials/Tutorial 6")
d <- read.csv("Data/kluane.csv")
d$plot <- factor(d$plot)
d$fenced <- factor(d$fenced,
levels = c(0, 1),
labels = c("Unfenced", "Fenced"))
d$fertilized <- factor(d$fertilized,
levels = c(0, 1),
labels = c("Unfertilized", "Fertilized"))
d$permanent <- factor(d$permanent,
levels = c(0, 1),
labels = c("Non-permanent", "Permanent"))
replications(phen.ach ~ fenced * fertilized * permanent + Error(plot), data = d)
fen <- as.fixed(d$fenced)
kluane.aov <- aov(phen.ach ~ fenced * fertilized * permanent +
Error(fenced:fertilized:plot), data = d)   # Fit the model
check_model(kluane.aov)
check_model(kluane.aov[[1]])
check_model(lm(kluane.aov))
library(GAD)
setwd("~/Library/CloudStorage/OneDrive-UniversityofNewBrunswick/Documents/Teaching/UNB Courses/2025-2026/Winter 2026/Experimental Design/Tutorials/Tutorial 6/Data")
install.packages(c("dbplyr", "dtplyr", "effects", "forecast", "ggiraph", "httr", "igraph", "later", "multcompView", "plotrix", "sp", "stars", "tseries"))
d <- read.csv("Data/kluane.csv")
d <- read.csv(kluane.csv")
d <- read.csv("kluane.csv")
d$plot <- factor(d$plot)
d$fenced <- factor(d$fenced,
levels = c(0, 1),
labels = c("Unfenced", "Fenced"))
d$fertilized <- factor(d$fertilized,
levels = c(0, 1),
labels = c("Unfertilized", "Fertilized"))
d$permanent <- factor(d$permanent,
levels = c(0, 1),
labels = c("Non-permanent", "Permanent"))
library(GAD)
fen <- as.fixed(d$fenced)
fer <- as.fixed(d$fertilized)
plo <- as.random(d$plot)
per <- as.fixed(d$perm)
kluane.lm <- lm(phen.ach~ fen * fer * per + # Main effects & interactions
(plo %in% fen:fer) + # Random plot effects
(plo %in% fen:fer):per, data = d) # Plot:permanence interaction
(kluane.tab <- gad(kluane.lm, quasi.f = TRUE))
comp.var(kluane.lm, anova.tab = kluane.tab)
d <- read.csv("kluane.csv")
op <- options(contrasts = c("contr.sum", "contr.poly")) kluane.aov <- aov(phen.ach~ fenced * fertilized * permanent +
op <- options(contrasts = c("contr.sum", "contr.poly"))
kluane.aov <- aov(phen.ach~ fenced * fertilized * permanent +
Error(fenced:fertilized:plot), data = d) summary(kluane.aov)
op <- options(contrasts = c("contr.sum", "contr.poly"))
kluane.aov <- aov(phen.ach~ fenced * fertilized * permanent +Error(fenced:fertilized:plot), data = d)
summary(kluane.aov)
d <- read.csv("Data/kluane.csv")
d <- read.csv("kluane.csv")
d$plot <- factor(d$plot)
d$fenced <- factor(d$fenced,
levels = c(0, 1),
labels = c("Unfenced", "Fenced"))
d$fertilized <- factor(d$fertilized,
levels = c(0, 1),
labels = c("Unfertilized", "Fertilized"))
d$permanent <- factor(d$permanent,
levels = c(0, 1),
labels = c("Non-permanent", "Permanent"))
op <- options(contrasts = c("contr.sum", "contr.poly"))
kluane.aov <- aov(phen.ach~ fenced * fertilized * permanent +Error(fenced:fertilized:plot), data = d)
summary(kluane.aov)
library(predictmeans)
residplot(kluane.aov)
pf(4.8841, df1=0.5172755*7, df2=0.5172755*133, lower.tail = F)
